/*
 * MatrixUtil.h
 *
 *  Created on: 2011-8-17
 *      Author: yong
 */

#ifndef MATRIXUTILE_H_
#define MATRIXUTILE_H_
#include "../Typedef.h"
using boost::numeric::ublas::matrix;
using boost::numeric::ublas::identity_matrix;
using boost::numeric::ublas::permutation_matrix;
namespace distrim
{
class MatrixUtil
{
private:
	typedef permutation_matrix<size_t> pmatrix;
	MatrixUtil() // Static utility class, no instance allowed.
	{
	}
public:
	/* Matrix inversion routine.
	 Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
	template<typename T>
	inline static bool LUInvertMatrix(const matrix<T>& input,
			matrix<T>& inverse)
	{
		BOOST_ASSERT(input.size1() == input.size2() && input.size1()
				== inverse.size1() && input.size2() == inverse.size2());
		// create a working copy of the input
		matrix<T> m(input);
		// create a permutation matrix for the LU-factorization
		pmatrix pm(m.size1());
		// perform LU-factorization
		int res = lu_factorize(m, pm);
		if (res != 0)
		{
			return false;
		}
		// create identity matrix of "inverse"
		inverse.assign(identity_matrix<T>(m.size1()));
		// backsubstitute to get the inverse
		lu_substitute(m, pm, inverse);
		return true;
	}

	/**
	 *@brief evaluate determinant of matix.
	 */
	template<typename T, typename R>
	inline static R Determinant(const matrix<T>& input)
	{
		BOOST_ASSERT(input.size1() == input.size2());
		matrix<T> m(input);
		pmatrix pm(m.size1());
		R det = 1.0;
		if (lu_factorize(m, pm) != 0)
		{
			det = 0.0;
		}
		else
		{
			for (size_t i = 0; i < m.size1(); i++)
			{
				det *= m(i, i); // multiply by elements on diagonal
			}
			det = det * detSign(pm);
		}
		return det;
	}

	template<typename T, typename R>
	inline static R Det(const matrix<T> &input)
	{
		Determinant(input);
	}

private:
	/**
	 * @brief determine sign of matrix determinat.
	 */
	static int detSign(const pmatrix& pm)
	{
		int sign = 1;
		size_t size = pm.size();
		for (size_t i = 0; i < size; ++i)
		{
			if (i != pm(i))
			{
				sign *= -1; // swap_rows would swap a pair of rows here, so we change sign
			}
		}
		return sign;
	}
};
}

#endif /* MATRIXUTILE_H_ */
